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We present an extensive pseudospectral study of the ran- 
domly forced Navier-Stokes equation (RFNSE) stirred by a 
stochastic force with zero mean and a variance ~ k 4 ~ d ~ v , 
where k is the wavevector and the dimension d = 3. We 
present the first evidence for multiscaling of velocity structure 
functions in this model for y > 4. We extract the multiscaling 
exponent ratios Cp/C2 by using extended self similarity (ESS), 
examine their dependence on y, and show that, if y = 4, they 
are in agreement with those obtained for the deterministically 
forced Navier-Stokes equation (3dNSE). We also show that 
well-defined vortex filaments, which appear clearly in studies 
of the 3dNSE, are absent in the RFNSE. 

PACS : 47.27.Gs, 47.27.Eq, 05.45.+b, 05.70.Jk 



Kolmogorov's classic work (K41) on homogeneous, 
isotropic fluid turbulence focussed on the scaling behav- 
ior of velocity v structure functions S p (r) = (|vj(x + r) — 
Vi(x)| p ), where the angular brackets denote an average 
over the statistical steady state He suggested that, 
for separations r = |r| in the inertial range, which is sub- 
stantial at large Reynolds numbers Re and lies between 
the forcing scale L and the dissipation scale rjd, these 
structure functions scale as S p ~ r^ p , with ( p — p/3. 
Subsequent experiments Q have suggested instead that 
multiscaling obtains with p/3 > £ p , which turns out to 
be a nonlinear, monotonically increasing function of p; 
this has also been borne out by numerical studies of the 
three-dimensional Navier-Stokes equation forced deter- 
ministically (3(fNSE) at large spatial scales |^||. The 
determination of the exponents £ p has been one of the 
central, but elusive, goals of the theory of turbulence. 
One of the promising starting points for such a theory 
is the randomly forced Navier-Stokes equation (RFNSE) 
||-||, driven by a Gaussian random force whose spatial 
Fourier transform f (k, i) has zero mean and a covariance 
< f i (k,t)f 7 -(k',t') >= Ak i - d -yp i] (\i)5{^ + \i l )8{t - f); 
here k, k' are wave numbers, t,t' times, i,j Carte- 
sian components in d dimensions, and -Py(k) the trans- 
verse projector which enforces the incompressiblity con- 
dition. One-loop renormalization-group (RG) studies of 
this RFNSE yield HJi| a K41 energy spectrum, namely, 
E(k) ~ k 2 S 2 (k) = fc 2 (|v(k)| 2 } - fc- 5 / 3 , if we set d = 3 
and y = 4; this has also been verified numerically ||. 
Nevertheless, these RG studies have been criticised for a 
variety of reasons |?],|^] such as using a large value for y in 
a small-?/ expansion and neglecting an infinity of marginal 
operators (if y = 4). These criticisms of the approxima- 



tions used in these studies might well be justified; but 
they clearly cannot be used to argue that the RFNSE 
is in itself inappropriate for a theory of turbulence. It 
is our purpose here to check if, indeed, the RFNSE is 
a good starting point for such a theory. Specifically we 
want to test whether structure functions in the RFNSE 
display the same multiscaling as in the 3<iNSE for some 
value of y; if they do, then we can argue that both equa- 
tions are in the same universality class and the RFNSE 
can, defensibly, be used to develop a statistical theory 
of inertial-range multiscaling in homogeneous, isotropic 
fluid turbulence. 

To achieve this end we have carried out an extensive 
pseudospectral study of the RFNSE and compared our 
results with earlier numerical studies |||| of the 3dNSE 
and experiments Q . We find several interesting and new 
results: We show that structure functions in the RFNSE 
display multiscaling for y > 4. As in the 3<iNSE, we ob- 
tain good estimates for ratios of the multiscaling expo- 
nents, such as Cp/C2j by using the cxtcndcd-self-similarity 
(ESS) procedure (Fig. 1) |Q; we obtain C 2 from S 2 (k) 
(Fig. 2). Next we investigate the y-dependence of ( p and 
find that it is close to the 3dNSE result (Fig. 3) for y = 4 
at least for p < 7. Thus the RFNSE should be a good 
starting point for a theory of inertial-range multiscaling 
in the 3<iNSE barring weak correction which do not af- 
fect ratios like Cp/Cz (see below). Furthermore we show 
that the qualitative behaviors of the probability distri- 
butions P(6v a (r)), where Sv a (r) = w Q (x) — w Q (x + r), 
are similar in the two models (Fig. 4). However, the 
shapes of constant- |w surfaces, where to is the vortic- 
ity, are markedly different (Fig. 5); the stochastic force 
destroys well-defined filamentary structures that obtain 
in 3<iNSE studies. This has implications for the She- 
Leveque fllf| formula for (p as we discus below. 

We use a pseudospectral method to solve the RFNSE 
numerically on a 64 3 grid with a cubic box of linear 
size L = 2ir and periodic boundary conditions; we have 
checked in representative cases that our results are un- 
changed if we use an 80 3 grid or aliasing. Aside from 
the stochastic forcing, our numerical scheme is the same 
as that used in Ref. |J. Our dissipation term, which 
is (v + i^ffc 2 )fc 2 v(k) in wave-vector (k) space includes 
both the viscosity v and the hyperviscosity vh; the ex- 
ponents are unaffected by vh if v > [pl]p|. Our 
numerical study of the RFNSE differs from conventional 
studies of the 3dNSE in two important ways: (I) For 
a fixed grid size we can attain higher Taylor-microscale 
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FIG. 1. Log-log plots (base 10) of S2(k) versus k for 
different values of y. The line indicates the K41 result 
<S2(fc) ~ fc _11//3 . The inset shows a repesentative plot of Re\ 
versus time(t) for y = 4. 

Reynolds numbers i?eA in the RFNSE, and hence a larger 
inertial range, than in the 3dNSE (Re\ ~ 120 compared 
to Re\ — 22 in our study), as noted earlier H for y = 4. 
(2) This advantage is reduced somewhat by the need to 
average statistical observables longer in the RFNSE than 
in the 3g?NSE. In the latter case it normally suffices to 
average over a few box-size eddy turnover times tl; this 
is not enough for the RFNSE since (a) Re\ fluctuates 
strongly over time scales considerably larger than tl (in- 
set of Fig. 2) and (b) the length of the f (k, t) time series 
required to obtain a specified variance for the stochastic 
force turns out to be quite large (the time series must 
be of length ~ 6tx to achieve the given variance within 
1 — 2%). Thus, in our studies, we have collected data 
for averages over 25 — 33tx (for different values of y ), 
after initial transients have been allowed to decay (over 
times ~ 10 — 20r^). Our T£ ~ 10r/, the integral-scale 
time used in some studies Jl2| ; tj = Lj/v rms , where the 
integral scale Lj = [J dkkE(k)/ J dkE(k)]^ 1 and v rms is 
the root-mean-square velocity. 

We begin by investigating the inertial-range scaling of 
the fc-space structure function 52 (fc) ~ k~^. Given this 
power-law form, the exponent £2 is easily seen to be re- 
lated to the r— space inertial-range exponent (2 by the 
formula (2 = ( 2 — 3. Our data in Fig. 1 for 4 < y 
are consistent with ( 2 = H/3 (i.e., the K41 value since 
E{k) ~ k 2 S 2 (k) ~ /c~ 5 / 3 ). For y = 4 this result has been 
reported earlier |(|. The y— independence of £ 2 above 
some critical y c (our data suggest y c ~ 4) is theoreti- 
cally satisfying since the variance of the stochastic force 
in the RFNSE rises rapidly at small k, so we might ex- 
pect that, for sufficiently large y, it approximates the 
conventional deterministic forcing of the 3c£NSE at large 
spatial scales. This point of view has been explored in 



the N — * 00 limit of an N— component generalization of 
the RFNSE 0; however, this study suggests C2 = 7 / 2 for 
y c = 4 < y; given our error bars (Table 1) it is difficult to 
distinguish this from the 0(y) RG prediction £ 2 = 11/3 
though our data are closer to the latter. For < y < 3 
both the one-loop RG || and the N — > 00 theory pre- 
dict £2 ~ l+2y/3+0(y ), which is in fair agreement with 
our numerical results, especially for small y (Table 1). We 
note in passing that, for < y < 4, there is no invariant 
energy cascade as in conventional K41: The dominance 
of dissipation at large k does lead to an energy cascade, 
but the energy flux depends on the length scale r; specif- 
ically II(r) w Ar v ~ 4 , with A the scale-independent part 
of the variance of the stochastic force in the RFNSE. 
A K41-type argument Jjl| now yields an energy-transfer 
rate ~ < Svf >/r ~ A y ~ A \ whence S^{r) ~ r( y ~ 3 ) and, 
if we assume simple scaling as in K41, Si,{r) ~ r^ -3 ) 2 / 3 , 
i.e., C 2 = 1 + 2y/3, which is same as the 0{y) RG pre- 
diction mentioned above. This formula breaks down for 
y < 0; however, the RG predicts correctly that the linear- 
hydrodynamics result obtains in this regime. 

Several precautions must be taken to ensure that sys- 
tematic errors do not affect the determination of C 2 - If 
k m ax is the largest wave- vector magnitude in our numer- 
ical scheme, we find that Lik max decreases with decreas- 
ing y; this shortens the inertial range of the energy spec- 
trum which can be used to obtain £ 2 . The lower the value 
of y the more difficult it is to obtain a dissipation range 
free of systematic, finite-resolution errors. For y < 4, we 
define kd = n^ 1 to be the inverse length scale at which the 
energy-transfer time t r ~ (r/v r ) ~ [Ar^ -6 )] 1 / 3 equals 
the diffusion time tn ~ [vk 2 + j/^/c 4 ] -1 ; this yields 

^ + Vh k d ^[Ak 6 d -y}^, (i) 

which when solved numerically shows that, for fixed A, 
kd increases as y decreases (Table 1). It is important to 
recognize that statistical steady states, with ill-resolved 
dissipation ranges that do not have a decaying tail [||, 
can be obtained by adjusting the amplitude A. In such 
cases kd ^> k max and we get spurious results for £ 2 . We 
find that, if we increase the hyperviscosity vjj, kd is suf- 
ficiently close to k m ax that we can resolve both inertial 
and dissipation ranges and obtain reliable values for £ 2 . 
Table 1 shows the range over which we fit our data for 
S2(k). Since our data for £ 2 indicate that y c ~ 4, we 
investigate multiscaling only for y > 4. 

Our data for £ 2 in Table 1 suggest that naive estimates 
for the multiscaling exponents £ p require longer inertial 
ranges than are available in our studies. However, we find 
that, as in the 3g?NSE, the extended-self-similarity (ESS) 
procedure (HQ 

can be used fruitfully here to extract 
the exponent ratios ( p /( q from the slopes of log-log plots 
of S p (r) versus S q (r) (see Fig. 2 for p = 5 and q = 2) 
since this extends the apparent inertial range. The ratios 
Cp/C2 that we obtain from such ESS plots are compared 
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FIG. 2. Log- log plots (base 10) of S$(r) versus <S2<V) illus- 
trating the ESS procedure. 
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FIG. 3. Inertial-range exponent ratios Cp/C2 versus p for 
the RFNSE with y = 4 and 6 (extracted from plots like Fig.2). 
The line indicates the SL formula. 

in Fig. 3 with the She-Leveque (SL) formula pJ| , which 
provides a convenient parametrization for the experimen- 
tal values for C, p . (If we assume the power-law form for 
52 (fc) in the inertial range (Fig. 1), we get Q 2 — C2 — 3 
and thence all the exponents Q p .) Figure 3 shows clearly 
that, with y = 4, our RFNSE exponent ratios lie very 
close to those for the 3dNSE and, to this extent, these 
two models are in the same universality class. The expo- 
nent ratios for y < 4 lie away from the 3<iNSE values as 
we might have anticipated from our results for 52 (k) (Fig. 
1). These results also lead to the expectation that, if y is 
sufficiently large, the exponent ratio Cp/C2 should become 
independent of y, for we find that £2(2/ = 4) ~ (2(2/ = 6). 
However, for p > 3, our data for Cp/Cziu = 6) fall system- 
atically below those for Cp/(2(y = 4) or the SL line. We 
also find that that the probabilty distributions of P(Sv r ) 
(Fig. 4) have non-Gaussian tails for r in the dissipation 
range; and for y > 4 the deviations from a Gaussian dis- 
tribution increase systematically with y. Thus, at least 
at the resolution of our calculation, it seems that the 
RFNSEs with y — 4 and y = 6 are in different universal- 
ity classes. However, we wish to point out that our data 
for y — 6 are more noisy than those for y — 4, so longer 
runs with finer grids might well be required to settle this 
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FIG. 4. Semilog plots of the distribution P(Sv r ) (i.e., 
P(5v a (r)) averaged over a) for r in the dissipation range and 
y — 4 and 6. A Gaussian distribution is shown for compari- 
son. 

issue conclusively. 

Strictly speaking the RFNSE with y = 4 falls in the 
same universality class as the 3<iNSE only in the ESS 
sense. For arbitrary y the energy flux through the k th 
shell is n fe = n(r = k- 1 ) ~ /* /£ (|f (k)| 2 )d 3 fc, where r 
is in the inertial-range and we have used Novikov's the- 
orem |H, i.e., (f(k) • v(-k)> - (|f(k)| 2 ). For y > 4, 
saturates to a constant for l;L » 1; but for y = 4, 
Tl k - log(fcL) in the RFNSE. This is to be con- 
trasted, with the 3(iNSE where = constant. Thus 
the inertial-range behaviors of all correlation functions 
in the two models are not the same. A K41-type dimen- 
sional analysis suggests that for y = 4 the energy flux 
Hit ~< 5v^ > jr ~ log(r/L); if we further assume that 
there is no multiscaling, then S p (r) ~ [r log(r/L)] p / 3 . 
Multiscaling will clearly modify this simple prediction; 
but some weak deviation from the Von-Karman-Howarth 
form 53 (r) ~ r must remain, since the standard deriva- 
tion of the Von-Karman-Howarth relation pl| does not 
go through @ with the RFNSE result for U k . To the 
extent that our data show that the ESS procedure works 
for the RFNSE, it seems that these weak deviations can- 
cel when we consider the ratios of structure functions; 
and, as noted above, for y — 4 the exponent ratios Cp/& 
agree with the SL result for the 3dNSE. 

Filamentary structures (Fig. 5) [0 in iso-|w| plots 
have been used as important ingredients in phenomeno- 
logical models for multiscaling in fluid turbulence. For 
example, the SL formula (ll) is obtained by postulat- 
ing a hierarchical relation among the moments of the 
scale-dependent energy dissipation; this yields a differ- 
ence equation for the exponents t p , which are simply 
related to the exponents C p ; one of the crucial bound- 
ary conditions used to solve this equation requires the 
codimension of the most intense structures. If these are 
taken to be vorticity filaments, their codimension is 2 
and one gets the SL formula. Filaments have been ob- 
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TABLE I. The dissipation-scale wavenumber kd (determined from Eq. 1), the integral-scale wavenumber ki = LJ 1 , the 
apparent inertial range over which we fit our data for S^fc), the hyperviscosities vh, the exponent £ 2 that we compute, and its 
0(y) RG value, for 1 < y < 4. The viscosity v is 5 x 10~ 4 in all these runs which use a 64 3 grid. 
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FIG. 5. Iso-|u>| surfaces obtained from instantaneous snap- 
shots of the vorticity fields showing filaments for the 3dNSE 
(left) and no filaments for the RFNSE with y = 4 (right). 



served in experiments also Jig]. We have shown above 
that the exponent ratios Cp/C2 that we obtain from the 
RFNSE with y = 4 agree with the SL formula. One might 
expect, therefore, that filamentary structures should ap- 
pear in iso-|u;| plots for the RFNSE. However, this is 
not the case as can be seen from the representative plot 
shown in Fig. 5. The stochastic forcing seems to destroy 
the well-defined filaments observed in the 3<iNSE without 
changing the multiscaling exponent ratios. Therefore, the 
existence of vorticity filaments is not crucial for obtain- 
ing these exponents, which is perhaps why simple shell 
models also yield good estimates for ( p . 

In summary, then, we have shown that the RFNSE 
with y = 4 exhibits the same multiscaling behavior as 
the 3dNSE, at least in the ESS sense. Probability distri- 
butions like P(Sv r ) (Fig. 4) are also qualitatively similar 
in the two models, in so far as they show deviations from 
Gaussian distributions for r in the dissipation range. It 
would be interesting to see if the RFNSE model can be 
obtained as an effective, inertial-range equation for fluid 
turbulence. We have tried to do this by a coarse-graining 
procedure that has been used p0| to map the Kuramoto- 
Sivashinsky(KS) equation onto the Kardar-Parisi-Zhang 
(KPZ) equation; however, it turns out that the 3dNSE 
— > RFNSE mapping, if it exists, is far more subtle than 
the KS — > KPZ mapping as we discuss elsewhere jn| . 
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